C     =============================================================70
C     Iteratively compute the heat transfer flux in DCMD process
      SUBROUTINE DCMD_FLUX(W, T, P, V, 
     &                    Re, Nu, TW, JH, JM, TPC, ThermEff, IFLAG)
C     Common variables declaration
      USE DCMD_CMOD
C     W - (input) tube-side and shell-side mass flowrates [kg/s]
C     T - (input) tube-side and shell-side temperatures [K]
C     P - (input) tube-side and shell-side pressure [Pa]
C     V - (input) tube-side and shell-side velocity [m/s]
C     Re - (output) Reynolds number
C     Nu - (output) Nusselt number
C     H - (output) film heat transfer coefficiency [W/m2-K]
C     JH - (output) heat transfer flux [J/m2]
C     JM - (output) mass transfer flux [kg/m2]
C     TPC - (output) temperature polarization coefficient
C     ThermEff - (output) thermal efficiency
C     IFLAG - (output) running status flag
      REAL*8 W(2), T(2), P(2), V(2),
     &       Re(2), Nu(2), TW(2), JH, JM, TPC, ThermEff
      INTEGER*4 IFLAG
C     Local variables declaration
C     rho - density [kg/m3], 
C     cp - heat capacity [J/kg-K], 
C     k - thermal conductivity [W/m-K],
C     mu - viscosity [Pa-s]
      REAL*8 rho(2), cp(2), k(2), mu(2)
C     d(1) - HF tube inner diameter [m], 
C     d(2) - HF tube outer diameter [m],
C     de(1) - tube-side equivalent diameter [m],
C     de(2) - shell-side equivalent diameter [m],
C     csa(1) - tube-side cross section area [m2],
C     csa(2) - shell-side cross section area [m2],
C     l - length of membrane [m],
C     c - permeability [kg/m2-s-Pa],
C     km - thermal conductivity [W/m2-K],
C     thick - thickness of membrane [m]
C     num - number of hollow fiber tubes
      REAL*8 d(2), de(2), csa(2), l, c, km, thick
      INTEGER num
C     local variables            
      REAL*8 G(2), Pr(2), Gr(2), VisCorr(2)
      REAL*8 TM, H(2), TWNEW(2)
      INTEGER i
C     function declaration
C     GradDP - correlate gradient of pressure drop [Pa]
C     GradVP - correlate gradient of water vapor pressure [Pa/K]
C     DHVAP - latent heat [J/kg]
C     SVP - saturated vapor pressure [Pa]
      REAL*8 GradDP, GradSP, DHVAP, SVP, SVP2
C     DRE - calculate Reynolds number
C     DPR - calculate Prandlt number
C     DGR - calculate Graetz number
C     DNU - calculate Nusselt number
      REAL*8 DRE, DPR, DGR, DNU

C     Retrieve common variables
      DO i = 1, 2
        d(i) = COM_ID(i)
        de(i) = COM_DEQ(i)
        csa(i) = COM_CSA(i)
        rho(i) = COM_RHO(i)
        mu(i) = COM_MU(i)
        cp(i) = COM_CP(i)
        k(i) = COM_LAMBDA(i)
      END DO
      l = COM_LEN
      c = COM_PERM
      km = COM_COND
      num = COM_NUM
      thick = COM_THK

C     Initialization
      IFLAG = 0

C     Computing
      DO i = 1, 2
C       Calculate flowing mass flux
        G(i) = rho(i)*V(i)
C       Calculate Reynolds number
        Re(i) = DRE(G(i), de(i), mu(i))
C       Calculate Prandtl number
        Pr(i) = DPR(cp(i), mu(i), k(i))
C       Calculate Graetz number
        Gr(i) = DGR(Re(i), Pr(i), de(i), l)
C       Calculate viscosity correction factor near the wall
        VisCorr(i) = 1.D0
C       Calculate Nusselt number
        Nu(i) = DNU(Gr(i), Re(i), Pr(i), VisCorr(i), i)       
C       Calculate film heat transfer coefficient
        H(i) = Nu(i)*k(i)/de(i)
      END DO
      
C     ***************************************************************
C     *                                                             *
C     *                    UNDER CONSTRUCTION                       *
C     *                                                             *
C     ***************************************************************      
      CALL DCMD_ESTW(T, H, TW, IFLAG)
      IF (IFLAG .EQ. 1) THEN
        TPC = (TW(1)-TW(2))/(T(1)-T(2))
        JM = c*(SVP(TW(1))-SVP(TW(2)))
        JH = JM*DHVAP(TW(1))+km/thick*(TW(1)-TW(2))
        ThermEff = JM*DHVAP(TW(1))/JH
      ELSE
        WRITE(*, *) 'Erron in wall temperature estimation'
        STOP
      END IF

!C     Iteratively compute the heat transfer flux
!      TW(1) = T(1)
!      TW(2) = T(2)
!C     Maximal number of iteration is set to 100
!      DO i = 1, 100
!        TM = 0.5D0*(TW(1)+TW(2))
!        JM = c*(SVP(TW(1))-SVP(TW(2)))
!        JH = JM*DHVAP(TM)+km/thick*(TW(1)-TW(2))
!        TWNEW(1) = T(1)-JH/H(1)
!        TWNEW(2) = JH/H(2)+T(2)
!        IF (DABS(TWNEW(1)-TW(1)) .LE. 1.D-5) THEN
!          IFLAG = i
!        ELSE
!          TW(1) = TWNEW(1)
!          IFLAG = -1
!        END IF
!        IF (DABS(TWNEW(2)-TW(2)) .LE. 1.D-5) THEN
!          IFLAG = i
!        ELSE
!          TW(2) = TWNEW(2)
!          IFLAG = -1
!        END IF
!        IF (IFLAG .NE. -1) THEN
!          IFLAG = i
!          TPC = (TW(1)-TW(2))/(T(1)-T(2))
!          ThermEff = JM*DHVAP(TW(1))/JH
!          RETURN
!        END IF
!        IF (i .GE. 100) THEN
!          WRITE(*, *) '! Error in wall temperatures searching'
!!          CALL WrtHis('! Error in wall temperatures searching',
!!     &                0, TWNEW, 0)
!          STOP
!        END IF  
!      END DO
      
      RETURN
      
      END SUBROUTINE